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With a view to developing a formalism that will be applicable at second perturbative order, we 
devise a new practical scheme for computing the gravitational self-force experienced by a point mass 
moving in a curved background spacetime. Our method works in the frequency domain and employs 
the effective-source approach, in which a distributional source for the retarded metric perturbation 
is replaced with an effective source for a certain regularized self-field. A key ingredient of the 
calculation is the analytic determination of an appropriate puncture field from which the effective 
source and regularized residual field can be calculated. In addition to its application in our effective- 
source method, we also show how this puncture field can be used to derive tensor-harmonic mode- 
sum regularization parameters that improve the efficiency of the traditional mode-sum procedure. 

To demonstrate the method, we calculate the first-order-in-the-mass-ratio self-force and redshift 
invariant for a point mass on a circular orbit in Schwarzschild spacetime. 


I. INTRODUCTION 

The detection of gravitational waves from inspiraling 
compact binary systems will be greatly assisted by theor¬ 
etical waveform templates. Using matched filtering tech¬ 
niques, these templates will help to extract the incoming 
signal from the detector noise and, once detection be¬ 
comes routine, the same templates will allow the para¬ 
meters of the source binary systems to be established. 
Constructing an appropriate bank of waveform templates 
necessitates understanding the two-body problem in the 
general relativistic context. Unlike its Newtonian coun¬ 
terpart, this problem cannot be solved analytically with 
closed-form solutions. Instead, a variety of methods ex¬ 
ist to approximate the solution, each best suited to a 
particular set of system parameters. 

To model binary systems of comparable mass and small 
orbital separation it is necessary to turn to numerical sim¬ 
ulations to solve the full non-linear Einstein field equa¬ 
tions. Numerical Relativity has made great progress in 
recent years and it is now routine to numerically evolve 
binary systems of comparable mass black holes or neut¬ 
ron stars through tens of orbits up to, and through, mer¬ 
ger [1, 2]. Aside from numerical truncation error, these 
results are exact. However, the computational cost of 
running simulations to merger increases rapidly as the 
mass ratio of the smaller to the more massive body is 
decreased. The computational cost similarly increases as 
the initial orbital separation increases. Consequently, in 
these domains other approaches are required. 

For systems with large orbital separation the post- 
Newtonian (PN) expansion can be employed. This is 
a perturbative approach to the problem that involves ex¬ 
panding the held equations in powers of the orbital ve¬ 
locity as a fraction of the speed of light. This approach 
has a long and rich history and today the post-Newtonian 


expansion of the dynamics of a binary is now known up 
to 3PN order, with many 4PN terms now known — see 
Ref. [3] for a recent review. 

When one of the bodies is substantially more massive 
than the other the small mass ratio of the system can 
be used as a perturbative parameter. In this approach 
the less massive of the two bodies is usually modeled as 
a point particle. Flux balance arguments allow the dis¬ 
sipative dynamics to be modeled [4, 5], but in order to 
include conservative corrections it is necessary to evalu¬ 
ate the local ‘self-force’ acting on the particle. With a 
point particle source this then necessitates a regulariza¬ 
tion procedure to remove the coulomb-like divergence of 
the metric perturbation that does not contribute to the 
orbital dynamics. Over the years this regularization pro¬ 
cedure has been placed on very firm theoretical footing at 
first order in the mass ratio [6-8] and recently has been 
understood at second perturbative order [9-11]. 

At first perturbative order a large number of practical 
self-force calculations have been made - see Ref. [12] for 
a recent review. Practical calculation techniques are of¬ 
ten prototyped with scalar-field toy models before the 
gravitational case is considered. In both cases motion 
in the spherically symmetric Schwarzschild spacetime is 
usually tackled first, before turning attention to motion 
in the more astrophysically relevant Kerr spacetime of a 
rotating black hole. 

In recent years, the calculation of gauge-invariant 
quantities has proven to be particular fruitful as it allows 
for comparison of self-force results with those of post- 
Newtonian theory and Numerical Relativity. A number 
of these gauge-invariant quantities are now known [13 — 
16] and using them to make cross cultural comparisons 
has been illuminating, helping to delineate the region in 
which perturbative approaches are valid, as well as act¬ 
ing as a cross-check on the widely different computational 
approaches taken by each scheme. One particularly in- 
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teresting result is that of Le Tiec et al. [17] showing that 
even for comparable mass systems perturbative meth¬ 
ods can make meaningful statements about the orbital 
dynamics. Calibrating Effective-One-Body theory is an¬ 
other important use for gauge-invariant results [18-20]. 

When making self-force calculations, an important 
consideration is the practical regularization technique to 
be employed. All the techniques derive from the same 
fundamental regularization procedure, but different tech¬ 
niques suit different calculational approaches. One of 
the most commonly employed techniques is the mode- 
sum prescription [21]. This approach relies on the mode- 
decomposition rendering the individual modes of the 
retarded field solution finite at the particle’s location 
and has had much success with 1+1 time-domain and 
frequency-domain calculations. For 2+1 or 3+1 time- 
domain decompositions, where the retarded held remains 
divergent at the particle, so called ‘effective source’ tech¬ 
niques were developed [22, 23]. This approach involves 
moving the contribution from the singular metric per¬ 
turbation near the particle into the source. This pro¬ 
cedure renders the otherwise divergent source finite and 
amenable to numerical treatment. 

To date all self-force calculations have been at first or¬ 
der in the mass ratio. The key motivation for the present 
work is to develop a set of techniques that can extend the 
work mentioned above to encompass second-order-in-the- 
mass-ratio calculations, with the aim of computing con¬ 
servative gauge-invariant quantities [24]. This immedi¬ 
ately suggests the basic form our approach should take: 

1. Work in the Lorenz gauge. Currently the regular¬ 
ization procedure at second order in the mass ratio 
is best understood in the Lorenz-gauge [9, 25]. 

2. Work in the frequency domain. At present it is not 
known how to stably evolve the monopole and di¬ 
pole contributions to the Lorenz-gauge linearized 
Einstein equation [26]. These instabilities are ob¬ 
served in 1 + 1D, 2 +ID and 3 +ID time-domain de¬ 
compositions on a Schwarzschilcl background, and 
similar instabilities have been seen in Kerr space- 
time [27]. Even if the multipole l > 2 modes are 
evolved separately in the time domain, the i = 0,1 
modes must be solved in the frequency domain. 

3. Regularize with an effective source. Unlike at first 
order in the mass ratio, the individual multipole 
modes of the second-order retarded field diverge at 
the particle’s location. This precludes the use of 
the mode-sum method for regularization as it re¬ 
quires the multipole modes of the retarded metric 
perturbation to be finite at particle’s location. 

The details of points 2 and 3 were fleshed out by the au¬ 
thors in Ref. [28] using a toy scalar-field example. In this 
work we address point 1 and extend our previous scalar- 
field results to cover Lorenz-gauge gravitational perturb¬ 
ations on a Schwarzschild background spacetime. 


In addition to laying out a formalism that will be ap¬ 
plicable at second order in the mass ratio, a natural by¬ 
product of this work is the extension of the standard 
mode-sum scheme to allow for direct regularization of 
the retarded tensor modes of the metric perturbation. In 
the standard mode-sum procedure the retarded tensor 
modes of the metric perturbation must be projected onto 
a basis of scalar spherical-harmonics before regularization 
can be performed. This step is cumbersome and, due to 
the coupling between scalar and tensor modes, requires 
the computation of additional tensor modes. Our new 
prescription neatly avoids these issues altogether. 

The layout of this article is as follows. In Sec. II 
we overview the Lorenz-gauge field equations and their 
retarded solution for circular orbits in a background 
Schwarzschild geometry. In Sec. Ill we discuss the stand¬ 
ard mode-sum and effective-source regularization proced¬ 
ures. In Sec. IV we detail the construction of an effective 
source for a Lorenz-gauge metric perturbation sourced 
by a particle on a circular orbit. In Sec. V we outline 
our numerical procedure and give our results. Finally, 
in Sec. VI we extend the mode-sum prescription to work 
directly with tensor-harmonic modes. There is additional 
supporting material in Appendices A-D. 

This paper follows the conventions of Misner, Thorne 
and Wheeler [29]; a “mostly positive” metric signature, 
(—,+,+,+), is used for the spacetime metric, the con¬ 
nection coefficients are defined by T x v = ^g Xa {g a ^,v + 
gav :l i - g^v,c r), the Riemann tensor is R T = F^ - 

T X^ + - r L r v> the Ricci tenSOr and SCalar 

are = R T tlT v and R = R,J J '. and the Einstein equa¬ 
tions are G^ v = R — ^g^R = 8 ttT^. Standard geo- 
metrized units are used, with c = G = 1. Greek in¬ 
dices are used for four-dimensional spacetime compon¬ 
ents and capital Latin letters are used for indices on the 
two-sphere. We work with standard Schwarzschild co¬ 
ordinates (i, r, 9 , cp) and also work with a second coordin¬ 
ate system, (t, r, a, /3), which is related by a rotation. We 
denote a point on the worldline of the point mass by a 
‘0’ subscript (e.g., cco) and indices on tensors evaluated 
on the worldline are indicated by an overbar (e.g. gp,p). 
We also find it useful to define / = 1 — 2 M/r where M 
is the mass of the background Schwarzschild black hole. 


II. LORENZ-GAUGE FIELD EQUATIONS AND 
THEIR RETARDED SOLUTION FOR A 
CIRCULAR ORBIT 

In this section we overview the Lorenz-gauge field 
equations for linear-in-the-mass-ratio perturbations of 
Schwarzschild spacetime, along with their decomposition 
into tensor-harmonic and frequency modes. The basis of 
tensor-harmonic modes we use is that of Barack and Sago 
[30, 31], themselves a modification of the basis given by 
Barack and Lousto [32] . Barack and Sago also further de¬ 
composed the monopole and dipole field equations into 
the frequency-domain to side-step instabilities that oc- 
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cur when evolving those modes in the time-domain (see 
Ref. [26] for a discussion of this issue). Later, the fre¬ 
quency decomposition was given for all modes with cal¬ 
culations for generic bound orbits by Akcay et al. [33, 34] 
and Osburn et al. [35]. 

A. Field Equations 

Let us denote by g the full spacetime metric, which we 
shall consider to be the sum of the metric perturbation, 
h, and the background Schwarzschild metric, g , such that 
g = g + h. Hereafter an over-ring will denote a quantity 
defined with respect to the background (vacuum) space- 
time. In a given coordinate system, the Einstein field 
equations will then take the form 

G^vig^v + hfn,] = 8-7T (2-1) 

where G is the Einstein tensor, a functional of the full 
spacetime metric <?, and T is the stress-energy tensor. Let 
us define the trace of the metric perturbation by Tr (h) = 
g^ v hfj, v . We shall find that the field equations for the 
metric perturbation take a simpler form when expressed 
in terms of the trace-reversed metric perturbation, h^ Vl 
defined by 

hfi V = Ti(h) (2.2) 

so named because Tr(/i) = — Tr(h). 

As discussed in the introduction, when we approach a 
second-order-in-the-mass-ratio calculation we will want 
to work in the Lorenz gauge. Consequently, to develop 
the necessary techniques we will work in the Lorenz- 
gauge with the first-order-in-the-mass-ratio calculation 
we present in this work. The Lorenz-gauge condition is 
defined by 

V^ = 0, (2.3) 

where the covariant derivative is taken with respect to the 
background metric. By expanding the Einstein tensor in 
powers of the mass-ratio and only retaining terms linear 
in /i we arrive at the (Lorenz-gauge) linearized Einstein 
equation given by 

Uh^ + 2 R p ° v h pa = -167 tT^ v (2.4) 

where □ = and R is the Riemann tensor of the 

background spacetime. In this work we shall take the 
metric perturbation to be sourced by a point particle of 
mass g. The corresponding energy-momentum tensor is 
given by 

/ GO 

[- det (g)]~ 1/2 S 4 (x fl - x%)u^u v dr , (2.5) 

-OO 

where det(g) = r 4 sin 2 0 is the determinant of the back¬ 
ground metric tensor, u M is the particle’s four-velocity 


and r is the proper time measured along the particle’s 
worldline. We also use to denote a general space- 
time coordinate and hereafter adopt then notation that 
a subscript ‘0’ denotes a quantity’s value evaluated at the 
particle. Note that for a circular orbit in the equatorial 
plane of a Schwarzschild black hole we have u r = ug = 0 
and u t = —E 0l = L 0 where 

Eo = fo \!w^3M' Lo=ro \lvo^3M' (2 ' 6) 

are the (specific) orbital energy and angular-momentum, 
respectively. Finally, we mention that the gauge equation 
(2.3) and field equation (2.4) are consistent so long as the 
particle is moving along a geodesic of the background 
spacetime (as then = 0). 

The field equation in the form of Eq. (2.4) is not well- 
suited to a numerical treatment as the metric perturb¬ 
ation diverges in the vicinity of the worldline. Instead, 
an effective-source approach can be employed to regular¬ 
ize the field equation and allow for a certain regular field 
to be solved for directly. Alternatively, with a 1 + ID 
or frequency-domain decomposition the individual mul¬ 
tipole modes of the metric perturbation become finite 
at the particle’s location and the mode-sum scheme can 
be employed to regularize on a mode-by-mode basis — 
see Sec. Ill below for an overview of these two regular¬ 
ization procedures. As discussed in the introduction, 
at second order it will become necessary to employ an 
effective-source scheme even within a multipole decom¬ 
position. For that reason, despite not being required for 
a first-order-in-the-mass-ratio calculation, we will pur¬ 
sue an effective-source approach within a multipole and 
Fourier decomposition of the field equations. We give the 
details of this decomposition now. 


B. Decomposition into tensor-harmonic and 
frequency modes 

In this section we overview the multipole and Four¬ 
ier decomposition of the metric perturbation and source. 
The explicit details of this decomposition have been 
laid out elsewhere [33, 36] and are summarized in Ap¬ 
pendix A; here we shall just present the key results re¬ 
quired for this work. 

There are many different conventions and notations 
used to define a tensor-harmonic basis. In this work we 
use the definition chosen by Barack and Lousto [36] with 
the slight modification introduced by Barack and Sago 
[30]. The key property of the Barack-Lousto-Sago tensor- 
harmonics is that they form a 10-dimensional basis for 
any second rank, symmetric 4-dimensional tensor field in 
Schwarzschild spacetime. This allows us to write the 10 
independent components of the (trace-reversed) metric 
perturbation in terms of the spherical-harmonic modes 
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of 10 fields, h^(t, r ) for * = 1,..., 10 via 

p2,7T n 7T 

^(*>r) = —CiT / / dtl (2.7) 

/id) Jo Jo 

where dfl = sin 0 dd dp and the details of the tensor basis 
(including definitions for and af*) are given in Ap¬ 

pendix A. We further decompose into Fourier-frequency 
modes, 

i r°° 

^m(*» r )= 2 ^y h { el(u,r)e- iut dt. (2.8) 


equations. This decomposition of the Lorenz-gauge con¬ 
dition is given by 



-/ 

(iu>mh^ + h ^ + —— 


(2.13) 

iuj m h {2) = 

- //#> + fhf 





r 

(h w - U 5) - 

a (3) - 

2A (6) ) , 

(2.14) 

iu m h {4) = 

_ / 
r 

[rUp + 2 h (5) 

+ Lh^ 

- *< 7 >) . 

(2.15) 

iu m h (8) = 

_ f 

r 

(rhW + 2 

- h < 10 >) 

5 

(2.16) 


For periodic motion the integral over frequencies reduces 
to a sum over discreet harmonics (hereafter modes). In 
particular, for a circular geodesic orbit the mode fre¬ 
quency is a simple overtone of the fundamental azimuthal 
frequency 


UJm — TnYly ,, 


(2.9) 


where = dp^/dt = \JM/r^ and m is the azimuthal 
mode index. In the remainder of this article, we can 
therefore denote the modes by ft^(r) without any ambi¬ 
guity. The stress-energy tensor can be similarly decom¬ 
posed into modes T^|(r) [32]. 

Substituting the mode expansions of and Xy into 
the linearized Einstein equation, Eq. (2.4), the angular 
and time dependence decouples. The spherical symmetry 
of the background geometry ensures that the individual 
multipole modes are eigenfunctions of the wave oper¬ 
ator (2.4) and consequently each multipole mode can be 
solved for independently from the others, though in gen¬ 
eral the ten tensorial components of each mode remain 
coupled. The resulting set of ordinary differential equa¬ 
tions for each multipole mode are given by 

n - 4/- 2 A4« (j) ^ = J&S(r - r 0 ), (2.10) 

where comes from the decomposition of the source 
and is the scalar wave operator, 

n l™ = ^2+y|;-r 2 [W-^]- (2.ii) 

Here, a prime denotes differentiation with respect to r 
and the potential term is given by 


Ve(r) = f{r) 


2 M 

r 3 


e(e + iy 

j.2 


( 2 . 12 ) 


The M (i \j) that appear in Eq. (2.10) are first order dif¬ 
ferential operators that couple the ten components of the 
metric perturbation; we give their explicit form in Ap¬ 
pendix B. In deriving the M. < ' l> ( :l j we use in this work, 
we have used the frequency-domain decomposition of the 
Lorenz-gauge condition (2.3) to simplify the resulting 


where L = £(£ + 1). 

The ten field equations (2.10) are not all coupled to¬ 
gether; instead they separate out into independent even 
(i = 1,..., 7) and odd (i = 8,9,10) parity sectors. Ex¬ 
amining the sources (given in Appendix B) we see that 

a dpt)]* =0 for £ + m = odd 

^0=8,9,10) K [yte( f / 2) Q v t)]* =0 for £ + m = even 

Consequently we have = 0 for £+m = odd and 

^ 0 = 8 , 9 , 10 ) _ o for l + rn = even. 

The gauge equations can be used to reduce the num¬ 
ber of fields that need to be solved for simultaneously. 
For example, for radiative modes (ui m ^ 0) in the odd 
sector one can solve for the for h^ and h ( 10 ) fields from 
which the h^ field can be constructed algebraically from 
gauge equation (2.16). Similarly, for radiative modes in 
the even sector the number of field equations to be solved 
simultaneously can be reduced by using Eqs. (2.13)- 

(2.15) . In this work we opt to only use Eqs. (2.14) and 

(2.15) to reduce the number of fields to be solved from 7 
to 5. The remaining gauge equation (2.13) can then be 
used as a consistency check on the final result. 

The static modes (ui m = 0) require a different treat¬ 
ment. In the odd sector both h^ and are zero as 
their sources vanish. The resulting equation for h ® can 
be solved for analytically - see Ref. [32] for details. In 
the even sector, gauge equations (2.14) and (2.15) can be 
used to eliminate the h ^ and h ^ fields which appear 
in Eqs. (Bl), (B3) and (B5). The resulting set of three 
ordinary differential equations were first solved numeric¬ 
ally [33], but more recently analytic solutions have been 
derived [35]. In this work opt for the numerical approach. 

For the non-radiative low multipole modes (£ = 
0,1, m = 0) analytic solutions are known [37]. The 
t = to = 1 mode is solved for numerically much as the 
other radiative even sector modes are, except for this 
mode h^ = 0 identically. We overview this hierarchical 
structure for solving the field equations in Table I. 


C. Retarded-field solution 

In this section we outline the calculation of the re¬ 
tarded field solution to Eq. (2.10) using the standard 
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l 

m 

£ + m = even 

l + m = odd 

0 

0 

(i) = l,3—>6 (A) 

- 

1 

0 

- 

(0=8 (A) 

1 

1 

(0 = 1,3,5,6-12,4 

- 

> 2 

0 

(0 = 1,3,5-A 6, 7 (A*) 

(0=8 (A) 

> 2 m ^ 0 

(0 = 1,3,5,6, 7 -A 2,4 

(*) = 9,10 -a 8 


Table I. Hierarchical structure for solving the field equations. 
A F implies the field(s) to the right should be algebraically 
constructed from the fields to the left using Eqs. (2.13)-(2.16). 
An (A) implies analytic solutions are known, and we employ 
them in this work except in the case of the even static modes. 


variation of parameters method. In this approach the in¬ 
homogeneous solution is constructed by multiplying the 
homogeneous solutions by suitable weighting coefficients. 
In each sector (odd/even, static/radiative) we must solve 
for k coupled equations and correspondingly the space of 
homogeneous solutions will be 2k dimensional. Using 
j = 1 ,,k as an index for the basis of homogeneous 
solutions, let us define the ‘inner’ and ‘outer’ homogen¬ 
eous solutions to the field equation by h^~ and h^ + , 
respectively. The inner solutions are regular at the hori¬ 
zon but diverge as r —> oo. Conversely, the outer solu¬ 
tions are regular at spatial infinity and diverge at the 
horizon. For the radiative modes the retarded solutions 
are selected by ensuring radiation at the horizon is purely 
in-going and radiation at spatial infinity is purely outgo¬ 
ing. This in turn implies that the asymptotic behaviour 
of the inner and outer solutions go as 


where the limits on the integral depend upon which 
weighting coefficient is being solved for. For the C ~’s 
a = r,b = o o and for the C'+’s a = 2 M, b = r. The source 
vector is given by k zeroes followed by the k sources from 
the right-hand side of the field equation (2.10). 

The delta function in the source means the integration 
can be done analytically and the inhomogeneous solu¬ 
tions can be written explicitly as 


/i (i) (r) 


where 


T, k j=iCpif + {r) r>r 0 
T, k j=i c Johf~{r) r<r 0 


Cl 

C, 


3 0 


4° = $ _1 M 


J (J) M 


( 2 . 21 ) 


( 2 . 22 ) 


Note that the C± 0 are r-independent constants. 


III. REGULARIZATION 

Building on the work of Mino, Sasaki and Tanaka [6] 
and Quinn and Wald [7], Detweiler and Whiting showed 
that the gravitational self-force can be computed as the 
derivative of a suitable regular metric perturbation, h^ v , 
via 


FUix o) = nk^ s Vsh* (x 0 ), (3.1) 


(r* —> ±oo) ~ e ± “ J ' r » ; (2.17) where 


where r* is the tortoise radial coordinate defined by 
dr „ I dr = /(r) _1 . A more in depth discussion of 
the asymptotic behavior of the radial fields is given in 
Refs. [33, 34], 

With the above definitions, the standard variation of 
parameters approach can be used to construct the in¬ 
homogeneous solutions to Eq. (2.10) via 


hM (r) = fa (r)~hf {r) + Cj{r)hf + {r) s j. (2.18) 
i=i 


To compute the weighting coefficients Cj we define a 
2 k x 2 k matrix of homogeneous solutions by 


$(r) 


1 

S3-* 

T 


i 

S-I 

T 

d r hf + 



(2.19) 


k^ s ^ u v v?u s 

+ ( 3 . 2 ) 

includes a projection operator that ensures the self-force 
is orthogonal to the particle’s four-velocity. 

The regular metric perturbation is constructed by sub¬ 
tracting an appropriate singular perturbation, h R IV , from 
the usual retarded metric perturbation h 1 ^, i.e., 

~KM = [h*A x ) - fyLO)] ■ (3-3) 

The construction of an appropriate singular field is dis¬ 
cussed at length in Refs. [8, 38]. One of the key features 
of the three metric perturbations is that they 

obey the field equations 


The weighting coefficients C~ (r) are then computed with 
the standard variation of parameters prescription: 



0 

J^\r’)5(r' - r 0 ) 


dr ', 

( 2 . 20 ) 


□ h^ s + 2 R\\hfJ S = —167rT^, (3.4) 

Dh^ / + 2R» l f v h% = 0, (3.5) 

from which we see that the retarded and singular perturb¬ 
ation diverge in the same way at the particle’s location, 
whilst their difference — the regular perturbation — is 











6 


smooth there. Using Eqs. (3.1) and (3.3) we can write 
the self-force as 


Keii( x o) = H lim V 4 (ft""( 

x—txo L 

= lim i F ?et( x ) ~ F s( X )b 

X^Xq 


(3.6) 


where 


self-force are given by 


C = E W - ^ [01 ) - ^ (3-10) 

£= o 

oo 

^eif = E (2i + 1) - - Fffi (3.11) 

1=0 


^ /s (*) = ^ 6 V s hiy s {x). (3.7) 

The divergence of at the particle makes it chal¬ 

lenging to work with Eq. (3.6) directly. Consequently, 
a number of reformulations have been devised to allow 
the gravitational self-force to be computed. Two of these 
schemes, the mode-sum method and the effective-source 
approach, we discuss now. 


A. Mode-sum method 


The key observation behind the mode-sum method is 
that although the full retarded and singular metric per¬ 
turbations are divergent at the particle, their individual 
multipole modes remain finite everywhere. The subtrac¬ 
tion between the retarded and singular contributions in 
Eq. (3.6) can then be made on a mode-by-mode basis. 
Explicitly we can write 


^eif(^o) = ^lim^ E [ F ret(aO - Ff(x) 
i=o 


(3.8) 


where a superscript £ denotes a quantity’s decomposition 
in to scalar spherical-harmonic modes and summed over 
m, i.e., 

„ £ f>‘2'K /*7T 

Kl ,S= Y l Y iJ*/ 2 ><P°)J 0 l F ret /8 YiJ0,V)dn. 
m——£ 

(3.9) 

We discuss below how we interface the tensor-mode com¬ 
putation of the retarded field outlined in Sec. IIC and the 
standard mode-sum scheme we are outlining now. 

The individual multipole modes of the retarded and 
singular contributions to the self-force, ir(ret/sjf are C~ l . 
That is, they are finite at the particle but, in general, 
their sided limits r —» Tq yield two different values, which 
we denote by F^, respectively. For circular orbits 
there is no closed-form analytic solution for Ff et , and 
typically it is computed numerically. The singular field, 
on the other hand, is amenable to an analytic treatment. 
The local structure of the singular field was first analyzed 
by Mino et al. [6] and Barack and Ori used these results 
to develop the mode-sum scheme shortly thereafter [21]. 

The scalar-harmonic mode-sum regularization formula 
for the redslrift invariant h^ u = h^u^u 1 ' [13] and for the 


The I-independent H^° 1 , , Ff, , D H , D M are known as 

regularization parameters and their value is known for 
generic geodesic orbits in Schwarzschild [39] and Kerr 
spacetime [40] . In general the coefficients of odd negative 
powers of l in the mode-sum formula are zero [41] and in 
the Lorenz-gauge D^ = Dh = 0. For circular orbits the 
other nonzero regularization parameters are given by 


ff[°] = 


Tr± 
[- 1 ] 


[ 0 ] 


4 A 4 

1C, 

(3.12) 

7 r\Ao + L o 



3M \ 1/2 
ro ) 

(3.13) 

Ao E o rc ot-i 

~ n(Ll + r 2 )V> [t ^ J ’ 

(3.14) 


where 1C = f^ 2 ( 1 — rQ M 2M sin 2 x ) 1 ^ 2 dx and £ = 

JJ r//2 ( 1 — sin 2 x) x l 2 dx are complete elliptic integ¬ 

rals of the first and second kind, respectively. 

The series in £ in both Eqs. (3.10) and (3.11) is trun¬ 
cated at £~ x . This is sufficient to regularize h uu and F r , 
but the resulting sum over £ converges rather slowly, with 
each term going as £~ 2 . It is possible to derive higher- 
order regularization parameters [41]; Ref. [42] provides 
the next two non-zero parameters that serve to increase 
the rate of convergence of the mode-sum to £ 6 . It is com¬ 
mon practice in mode-sum calculations to numerically fit 
for the yet higher-order unknown parameters to further 
increase the rate of convergence of the mode-sum. 

Lastly, as we mentioned above, we compute the re¬ 
tarded metric perturbation within a tensor-harmonic de¬ 
composition, whereas the standard mode-sum approach 
requires the retarded metric perturbation decomposed 
into scalar-harmonic modes as input. Thus before reg¬ 
ularizing we must project the tensor-harmonic modes of 
the metric perturbation onto a basis of scalar harmonics. 
The projection equation takes the form 


-^ret 




i 3 


£ £ yH*/ 


(3.15) 


where the details of the (but not the self-force 

obtained after summing over £) depends on the way in 
which the definition for the force (a quantity which is 
defined on the worldline) is extended off the worldline 
to the whole two-sphere. Barack and Sago [30] made the 
computationally-convenient choice of k abcd being given 
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by u a having a constant value on the two-sphere, and the 
metric having its usual tensorial value. Their expressions 

for the are rather cumbersome so we do not give 

them here; instead, their explicit form can be found in 
Appendix C of Ref. [30] . Likewise, a similar formula can 
be derived for h uu [43]. 

The sum over p in Eq. (3.15) means that in order 
to compute the self-force by regularizing £ max scalar- 
harmonic modes one must compute (£ max + 3) tensor- 
harmonic modes. Similarly for h uu one must compute 
(f max + 2) tensor modes. In Sec. VI below we will re¬ 
cast the standard mode-sum formula to use tensor modes 
rather the scalar modes, which will avoid this projection 
step altogether. 


B. Effective-source approach 

The effective-source approach is an alternative prac¬ 
tical regularization scheme for handling the divergence 
of the retarded field. Rather than first computing the 
retarded field and then subtracting the singular piece as 
a post-processing step, as in the mode-sum scheme, one 
can instead work directly with an equation for the reg¬ 
ular field. This idea was first proposed in Refs. [22, 23] 
and has the distinct advantage of involving only regu¬ 
lar quantities, making it applicable in a wider variety of 
scenarios than the mode-sum scheme. 

Using Eq. (3.3) to rewrite h T °* in terms of h^ v and , 
we can rewrite Eq. (2.4) as 

+ 2 R\\h R pn = -16t tT^ - Dh% - 2 R p ° v h% 

(3.16) 

If h% is precisely the Detweiler-Whiting singular met¬ 
ric perturbation, then the two terms on the right hand 
side of this equation cancel and h^ v becomes a homo¬ 
geneous solution of the wave equation. However, one 
typically does not have access to an exact expression for 
hp V as the Detweiler-Whiting singular metric perturba¬ 
tion is defined through a Hadamard parametrix which is 
not even globally defined. Instead, the best one can typ¬ 
ically do is a local expansion which is valid only in the 
vicinity of the worldline. Let us denote such an approx¬ 
imation to h^ v by h^ v . With the latter we will construct 
an effective source that will allow us to directly compute 
the regular field at the worldline. 

The puncture field h^ v is only valid near the worldline 
and so, to avoid ambiguities in the definition of the effect¬ 
ive source, one must ensure that the puncture field goes 
to zero far from the particle. This is can be achieved by 
multiplying h^ v by a window function, W, with proper¬ 
ties such that multiplying it by h’f only modifies terms 
of higher order in the local expansion about the worldline 
than those which are explicitly given in h^. In our par¬ 
ticular case, it suffices to choose W such that W(x 0 ) = 1, 
W(xq) = 0, W"(:to) = 0 and W = 0 far away from the 


worldline. The residual metric perturbation, hffjf, then 
obeys 




2 R p °«h% = S eS , 


(3.17) 


where the effective source is given by 

S eff = -16 nT^ - □(W^„) - 2 Rr; v (yVh p p<y ). (3.18) 

This effective source is smooth and finite everywhere, but 
has limited differentiability on the worldline. The corres¬ 
ponding residual field has the properties 

= h* v (x 0 ), V s h™(x 0 ) = V s h* v (x o), 
h™{x) = h™*(x) for x £ supp(W). (3.19) 

As the residual metric perturbation coincides with the 
retarded metric perturbation far from the particle we 
can use the usual retarded metric perturbation bound¬ 
ary conditions when solving Eq. (3.17). 


IV. EFFECTIVE SOURCE IN THE 
FREQUENCY DOMAIN 

A. Construction of the puncture fields 

1. Coordinate expansion of the singular field 

At the core of our calculation is an effective source 
for the field equations which is constructed from an ap¬ 
proximation to the Detweiler-Whiting singular field. A 
suitable covariant expansion of the singular field is given 

by 


h ab = 9a a gb 


1 UaUl 


0(e) 


(4.1) 


where e is an order-counting parameter, s = (g a i + 
Ua,ui)<J a <J b , u a is the four-velocity and g a i the background 
metric, with both defined as tensors on the worldline (i.e. 
at the spacetime point x 0 ). We have also introduced the 
bi-vector of parallel transport g a a {x,x o) and the Synge 
world-function cr(x,Xo), both of which are functions of 
the worldline point Xq and the point where the singular 
field is to be evaluated, x. 

This approximation is sufficient to produce a residual 
field which is finite on the worldline and which gives 
the correct, regularized self-force. Several higher order 
terms in this expansion are also known [42] and can 
be incorporated into the calculation in order to accel¬ 
erate convergence. However, for clarity we illustrate 
the approach with this simple low-order approximation 
and note that the methodology does not fundamentally 
change at higher orders. 

We now wish to use the approximation (4.1) as a start- 
ing point to compute the puncture fields h)f n . To this 





end, we follow previous regularization strategies by in¬ 
troducing a Riemann normal coordinate system in the 
vicinity of the worldline, and rewrite (4.1) as a coordin¬ 
ate expansion in terms of these coordinates. Specific¬ 
ally, we assume that the spacetime can be represented in 
terms of a spherical coordinate system with polar and azi¬ 
muthal coordinates a and /?, radius r and time t. Note 
that, although our focus here is on the Schwarzschild 
spacetime, the assumption of a spherical coordinate sys¬ 
tem does not necessarily limit us to spherical symmetry; 
for example, the method works equally well in the non- 
spherically symmetric Kerr spacetime [44]. 

Now, orienting our coordinate system such that the 
worldline is instantaneously at a = 0, we define the 
Riemann normal coordinates w\ = 2 sin j cos /3 and W 2 = 
2 sin | sin /3. Using coordinate expansions of g a a {x,x 0 ) 
and a(x,Xo) about x = xq to linear order in x — Xq, 
we obtain an approximation to (4.1) in terms of the 
(t,r,wi,W 2 ) Riemann normal coordinate system. Struc¬ 
turally, our coordinate expansion has the form 


r(S) _ 1 C ab 

hab - e P 


C ab Ar 


c^Ar 3 ! 


O(e), (4.2) 


where p is the leading-order term in the coordinate ex¬ 
pansion of s and the coefficients and do not 

depend on A r or a (and hence w\ and W 2) 1 ■ The coef¬ 
ficients are also independent of t since we have chosen 
At = 0, i.e., x and Xq are points on the same time 
slice. There is still a potential time dependence, how¬ 
ever, through the dependence of the coefficients on the 
worldline and four-velocity. 

In the next subsection, we will seek a decomposi¬ 
tion into spherical-harmonic modes. We therefore apply 
the (approximate) Jacobian from (wi, W 2 ) coordinates to 
(a,/3) coordinates. In doing so, we pull out a factor of 
sin a from the Jacobian when computing h t p, h r p and 
h a f 3 , and a factor of sin 2 a when computing hpp. The 
reason for doing so will become clear during the mode de¬ 
composition, and is related to the fact that the Riemann 
normal coordinate system is regular on the worldline, but 
the (a,/3) coordinate system is not. 

Evaluating Eq. (4.1) for our particular case of a circular 
orbit in Schwarzschild spacetime, we arrive at our desired 
coordinate expansion of the Detweiler-Whiting singular 
metric perturbation. With Riemann normal components 


1 This form is valid for the case of circular orbits in Schwarzschild 
spacetime, where any quadratic dependence on w 1 and W 2 can 
be replaced with a term involving p 2 and A r 2 . The structure 
is slightly more complicated in more general cases where odd 
powers of wi and W 2 can appear, but nonetheless the following 
analysis remains qualitatively unchanged. 


given by 


htW! 


4roU v ,(ro — 2AT) 2ArroU ¥ 


pi ro — 3AT ro — 3 AT 

— 3Afro + 2 AT 2 — 2AT 2 sin 2 /3 


forwi — 


(r 0 -2Af)(l-^| 37 sin 2 /3) 
4ATrosina cos/3 
p( r o ~ 3AT) 
cos 2 /3 


4 Afrp 


2Ar Afro 
p Lro — 3AT ' ?'o — 3Af 
3r 0 - 7AT - 2AT sin 2 /3 


(r 0 - 2Af)(l - Abj siir /?) 


(4.3a) 

(4.3b) 

(4.3c) 


our approximation to the Detweiler-Whiting singular 
metric is then given by 


htt — — 


4(r 0 - 2Af) 2 


2Ar 


P Lr 0 (ro — 3AT) r5(r 0 - 3AT) 
r% - 7Afr 0 + lOAf 2 - 2Af (r 0 - 4AT) sin 2 /3 


1 - 


M 


r 0 -2 M 


sin /3 


ht r 


4roD ¥ ,(r 0 — 2Af) sinct cos/3 


p(r 0 - 3AT) 
hta COS /3, 

ht /3 = -h tWl sin a sin /3, 

— 0, 

fora — forwi COS/3, 
fo r p = —h rw r sin a sin /3, 
foaa — fow\w\ COS /3, 
hap hit 


sin a sin /3 cos /3, 


him = hr„ 7 ,), sin 2 a sin 2 /3. 


(4.4a) 

(4.4b) 

(4.4c) 

(4.4d) 

(4.4e) 

(4.4f) 

(4.4g) 

(4.4h) 

(4.4i) 

(4-4j) 


This approximation includes all contributions at order 
e^ 1 and e°, with the exception of terms proportional to 
Ar 3 /p 3 , which we neglect as their mode decomposition 
yields only terms proportional to Ar 2 and higher. 


2. Mode decomposition 


We now proceed with the decomposition of our co¬ 
ordinate expansion into tensor spherical harmonic modes. 
For this, we must evaluate the integrals of the singular 
field against the tensor spherical harmonics, 


jsv p - r 

n em ~ 


pa} 


W 


1 TK<\ '/ 1 


‘ sin a da d/3. 


(4.5) 

For the circular orbit case we are considering here, the 
explicit form for the integrand for each i = 1 ,..., 10 field 
is given in Table II. 

The mode decomposition works much the same as with 
the scalar field case described in [28]. There are, however, 
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Integrand 


(M,0) 
(M, ±1) 
( 3 , 40 ) 
( 4,4 ± 1 ) 


r \j^r{htt + f 2 hrr)Pj( cosa) sin a 

±2r \J ii$i^Tjf htr cos P p e( cos a ) sina 

r^J^^fjhtt - f 2 h rr )P°(cos a) sin a 
±2 sj (2t+1 \^ t+1) [ sin 2 /3 Pi {cos a) + cos 2 4 P<+l( 


2£+l 


(5,40) 


2 \/~4irfhrini cos^^r [-Pj+i^osa) - P°i(cosa)j 


(5,4±2) 2y/ cos ^W+T) [4sin 2 fiP 2 (cosa) + (cos 2 /3 - sin 2 /3) 


(Z—l)ZP£ +1 (cos at) —(Z+1)(Z+2)P£_ 1 (cos o:) I 
2£+l 


(6,40) 


-\J ^rh wiwi p°( cos a) sin a 


(7,4 ±2) 


4 


2 . /(2l+l)(l-l)l(l+l)q+2) 


47T 


(£-l)£(£+l)(£+2) sin 

\2 n2 


[ 8 COS 2 /3 Sin 2 /3 ( f - 1 ) 2p i+l( c ”°0-d+2) 2 P/U(c OSa ) 


1 C„„„2 a „ ; „2 ^2 (^—1) (2£— l)P^_ 2 ( cos a)—2(£—3)(£— 1 )(r+2)(^+ 4 )(2r+ 1 )P^(cos a) + (r+!) 2 (r+2) 2 (2r+3)P^_ 2 (cos a) / 

+ (COS p Sin P) 2(21—1) (21+1) (21+3) | 


( 8 , 4 ± 1 ) - 2 \] {2i+1 li (i+1) - w^) [^(cosa) cos 2 0 


Z^ P}, i (cos at) — (£+l)^P/_i (cos a) 

+ - ± - 2Z+1 - sm 


m 2 /?] 


(9,4 ±2) T \j Al{t-l) lJ^) Ai f h ^i cos Pl(jtpr) [( cos2 P - sin2 0)Pe(cosa) + sin 2 (5 (t 1)£Pg + 


(cos a:) — (£+l)(£+2)Pj^ ^ (cos a) 1 


2Z+1 


( 10,4 ± 2 ) 


,4^ 


)(l-l)l(l+l)(l+2) 


lih .‘ l;iu ’ 1 [(cos 2 /3 — sin 2 /3) 

sin a v ' 


2 (£— l) 2 P4_i(cos a) —(£+2) 2 P 2 _ 1 (cos a) 


(e-i)e(e+i)(e+2) 

2 o ■ 2 0 {t-\) 2 1 2 (2t-\)pj, 2 (cosa)-2(r-3)(i-l)(r+2)(<+4)(2r+l)P?(co S a) + (r+l) 2 (r+2) 2 (2^+3)P?_ 2 (co Sa 

+ cos p sm p -^- 


(2£-l)(2£+l)(2£+3) 


(I 


Table II. Integrands appearing in the mode decomposition of all 10 tensor-harmonic components of the singular metric 
perturbation for the case of a circular geodesic orbit in Schwarzschild spacetime. 


some key differences which introduce additional complex¬ 
ity to the gravitational case: 

1. The fact that we have tensor (as opposed to scalar) 
harmonics makes the mode decomposition integrals 
slightly more involved. 

2. Whereas in the scalar case we only required the 
m! = 0 modes, we now require m! = 0 for 44 4„ 
and 44 m' = 1 for 44 4m and 44 m! = 0,2 
for 44 and to' = 2 for h £>, 4m and 44- This is 
because we would like to compute the metric per¬ 
turbation and its derivative (for the self-force) on 
the worldline, and those are the only modes which 
don’t vanish on the worldline, at a = 0. Note that 
in principle other modes could contribute (to' = 1 

for 44 and 44 m' = 0 for 44 to' = 0,2 
for 4m and 44 to' = ^1 for 42 and 42 > and 

to' = 1,3 for 42 and 4m 2, but the integrals for 
those modes all contain odd powers of sin /3 or cos /? 
and therefore their contribution vanishes after in¬ 


tegration over /3. 

3. The coordinate approximation we are using for the 
singular field has a spurious non-smoothness away 
from the worldline, at a = ir. This can be seen 
in p , which has a /3-direction dependent limit as 
a —> ir. This problem did not manifest itself in the 
scalar case, since the isotropic nature of the m' = 
0 mode means it cannot include any information 
about direction dependence. 

The first two items above do not cause any funda¬ 
mental issues, they merely add some extra algebraic 
complexity to the problem. The third item, however, 
does cause problems if not handled appropriately. The 
non-smoothness introduces a spurious component in the 
puncture which behaves as v ,' in a mode-sum formula 
such as Eqs. (3.10) and (3.11). This renders the sum not 
absolutely convergent, although the (—1) £ factor means 
that it is in fact conditionally convergent since, for ex¬ 
ample, (~4(4i) f+1 ) = 1- In practice, this makes 

the sum over modes converge very slowly, see Fig. 1. 
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Conditional and absolute convergence of the mode-sum 
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Figure 1. Effect of a spurious non-smoothness away from 
the worldline on the convergence of a mode-sum scheme near 
the worldline. The ^non-smoothness manifests itself as a term 
of the form , and appears to spoil any hope of 

rapid convergence (blue/orange dots). This can be mitigated 
by using a smoothing factor which converts the conditionally 
convergent behaviour into a more rapid absolutely convergent 
behaviour (green). The final result is not altered since the 
infinite sum of conditionally convergent terms is exactly equal 
to the infinite sum of absolutely convergent terms. 


Fortunately, there is a straightforward resolution to 
this problem. A smooth window function, W m ' (a), in the 
a direction is effective in eliminating the spurious non¬ 
smoothness affecting the modes. To ensure the self-force 
is not affected, we require that W m ' («) ~ 1 + a 2 near 
a = 0, while eliminating the effect of the non-smoothness 
on a particular m ' mode requires W m '(a ) ~ ( 7 r —a)T m / 2 1 
near a = 7 r, where \m! /2] is the smallest integer greater 
than or equal to m' /2. We make the particular choice 
W m > (a) = (cos §) I " m , which satisfies both of the 
above criteria. 


3. Integrals over a 

In the circular geodesic case, the quantity p appearing 
in the singular metric perturbation is given by 


P = 


2x r o( r o - 2M) 2 


vq — 3 M 


(5 2 + 1 — cos a), 


where 


and 


6 2 = 


A?’ 2 ro — 3 M 

2xro(ro-2M)^ 


x= 1- 


M 


ro — 2 M 


sin 2 /3. 


(4.6) 


(4.7) 


(4.8) 


Then, the integrals over a all take one of nine possible 
forms which can be evaluated analytically. In our par¬ 
ticular case, we are only interested in the behaviour at 


the leading two orders in Ar = r — r 0 . To simplify our 
expressions, we introduce 


A — A _ 
Ai = Ai ; q — 


1(1+1) 


and 


A 2 = A 2 .fi = 


(2£-l)(2£+3) > 


(£—l)£(£+l)(£+2) 

2.0 (2£-3)(2£-l)(2£+3)(2£+5) • 


(4.9) 


(4.10) 


Then, for and the Ar-expanded integrals 

are given by 

Pj, } (cos a) sin a 

Jo (<5 2 + 1 — cos a ) 1 / 2 


da 


1 


21+11 


For hff, they are given by 

P\ (cos a) sin 2 a 


2V2 ~ 2{2l + l)\5\ +e>( 6 2 )l. (4.11) 


/0 (S 2 + 1 — cos a ) 1 / 2 
1 


da 


21+ 1 


— 8a/2Ai + 0(S 2 ) , (4.12) 


For and they are given by 

r Wi(a) t 2 Pl +1 {cosa) - (£ + l) 2 P/_ 1 (cosa) 
Jo £(£ + 1 ) (2£ + 1)(<5 2 + 1 — cos a ) 1 / 2 

1 r 6x^2 


da 


2£ + 1 L (2£-l)(2£ + 3) 


+ (2I+l)\8\+0(5 2 


(4.13) 


and 


H’i(a) P }{cos a) 


da 


/0 £(£ + 1) (<5 2 + 1 — cosa ) 1 / 2 

1 r . 6\/2 


2i+ 1 L 


- 8v^Ai + 


{21 — 1)(2£ + 3) 
+(2£+l)|<5| + 0(A 2 )1. (4.14) 


For h\ 0 they are given by 

['*£(£ + 1) [P ° +1 (cos a) - P°_ 1 (cosa)] sin a 


/0 (2*+l) 
1 


(8 2 + 1 — cos a) 1 / 2 

-8v^Ai + 0(8 2 ) . 


21 + 1 L 

For h ^2 and they are given by 

/’ 7r 1 W 2 (a)P 2 {cosa) sina 

Jo £{£ + 1) (8 2 + 1 — cos a) 1 / 2 

= —-— [ 32 V 2 A 2 
2 £ + 1 L 


da 

(4.15) 


da 


V2Qy/2(t-!)(£ +2) 

' (2 1 - 3)(2£ - 1)(2£ + 3)(2£ + 5) 


0(8 2 ) \, 

(4.16) 
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and 


W 2 (a:) sin a 


/o (2£ + 1)£(£ + I)(i5 2 + 1 — cos a) 1 / 2 

(£ — l)FP 2 +1 (cosa) — (£ + 1 )(£ + 2)P 2 _ 1 (cosa) 
1 


da 


2 £ + 1 L 

+ 


- 32y/2A 2 
240v / 2(f — !){£ + 2) 


(21 - 3)(2£ - 1)(2£ + 3)(2£ + 5) 


0(6 2 ) .(4.17) 


Finally, for and h they are given by 

I"* W 2 (a) esc a 

Jo (£ - !)£(£ + !)(* + 2)(2£ + l)(d 2 + 1 - cos a) 1 / 2 


£— l)“P/ +1 (cosa:) — (£+ 2) 2 P/_ 1 (cosa) da 


(4.18) 


r() M 0M ■ All other cases can be related to these three using 
the recursion relation for the hypergeometric function, 


T p+ \(k) — _ ^.Fp-i(fc) + 

where F p (k) = 2 Pi (p,^,l,k). 


1 -2p+ (p- \) 

P {k - 1 ) 


-Pp(fc), 

(4.21) 


B. Construction of the effective source and 
residual fields 

In order to construct an effective source we must choose 
a window function, W, to confine the definition of the 
puncture to a neighbourhood of the worldline. As dis¬ 
cussed in Sec. Ill B the constraints on the window func¬ 
tion are that W(xo) = 1, >V'(a; 0 ) = 0, W'ixo) = 0 and 
W = 0 far away from the worldline. These conditions 
leave considerable freedom when choosing a window func¬ 
tion. In this work we shall use the window function given 
by 


and 


/o (£-l)£(£ + l)(£ + 2) 
W 2 (a) esc a 


(6 2 + 1 - cos a) 1 / 2 (2£ - l)(2£ + l)(2£ + 3) 

(£ - 1) 2 £ 2 {2£ - l)P/ +2 (cosa) 

—2(£ — 3)(£— 1)(£ + 2)(£ + 4)(2£ + l)P { 2 (cosa ) 


+(£ + lf(£ + 2) 2 (2£ + 3)P 2 _ 2 (cos a) 
1 r nz . 40^2 


da 


2£ + ll 


- 32v / 2A 2 + 


(2t-l){2£ + 3) 

+(2£+l)\5\+0(5 2 )\. (4.19) 


4 . Integrals over /3 

With the integrals over a having been evaluated ana¬ 
lytically as a power series in 6, we are next faced with 
the integrals over /3. The functional dependence on /3 can 
be rewritten in terms of integer and half-integer powers 
of x = 1 — ro 0 2M sin 2 jd. These integrals are straightfor¬ 
ward to evaluate and yield either polynomials in ro M 2 M > 
or complete elliptic integrals with argument ro M 2M • Spe¬ 
cifically, 


/ x " d (d = 2n 2 F 1 (n,ll l7 ^ M ) 1 (4.20) 

J 0 

which has three special cases: for n = —1/2 it reduces 
to the elliptic integral of the first kind, /C( ro /^ M ); for 
n = 1/2 it reduces to the elliptic integral of the second 
kind, £( ro /^ M ); for n an integer it is a polynomial in 


W(r) = e~ 8M 4(r-ro)4 . (4.22) 

We make this choice as it is easy to implement and, al¬ 
though not formally compact, it is effectively compact 
within our numerical scheme. Other authors have made 
different choices. Vega et al. [45] used a compact window 
function that allowed for a smooth transition from the 
residual to the retarded held. Alternatively, a compact 
source can be achieved using the worldtube approach of 
Barack and Golbourn [22]. In Ref. [28] we used a Heav¬ 
iside n function and showed that this was equivalent to 
the worldtube method. In this work we opt not to do this 
for ease of implementation, though we note, by building 
on a draft of this work, that a worldtube method has 
been implemented for the gravitational case [46]. 

With the window function chosen the effective sources 
are given by 

S^ S =J&8(r - r 0 ) - (w~hf m P ) 

+ 4 r 2 M^ {j) {wh^l P ). (4.23) 

For brevity we will not display the explicit form of the 
5« eff Using the Held equations in Appendix B and 
punctures in Appendix C it is straightforward to com¬ 
pute the effective sources using computer algebra pack¬ 
ages. However, we do point out one potential subtlety: in 
the above equation we have implicitly assumed that the 
wave operator commutes with the mode decomposition, a 
fact which is not necessarily true. Indeed Barack and Ori 
[47] pointed out that the mode decomposition does not 
always commute with radial derivatives; likewise from 
the Wigner-Eckart theorem one may be concerned that 
a spherical-harmonic mode decomposition which fails to 
include all modes would not commute with the angular 
derivatives. In the current context both concerns turn 
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out to be unfounded. The Barack-Ori observation is only 
an issue if the limit A?’ —► 0 is taken, but we avoid doing 
so while computing the puncture fields. The neglection 
of higher spherical-harmonic modes of the puncture does 
indeed affect the effective source we obtain, but only in 
a way which affects the higher derivatives of the residual 
field (since those higher modes vanish when evaluated at 
a = 0). 

The construction of the residual metric perturbation 
now proceeds as follows. Via the variation of parameters 
prescription we have 

k 

h^ ies (r) = (o- Tes (r)hf-(r) + Cf ies (r)hf + (r)') , 

j=i 

(4.24) 

where recall we use j to index the k basis of a given £m 
mode. The weighting coefficients are given by 

( cf^r) )=£ V) ( ) *•'> ( 425 ) 

where $ is the 2k x 2k matrix of homogeneous solutions, 
defined in Eq. (2.19). The source vector is formed of k 
zeroes followed by the k effective sources. The integration 
limits in Eq. (4.25) depend upon which weighting coeffi¬ 
cient is being solved for. For the C“ res ’s a = r, b = oo 
and for the C') l ~ res ’s a = 2M, b = r. 

In order to compute the self-force we also require the 
first radial derivatives of the metric perturbation fields. 
These are easily constructed via 

k 

h^ ies '{r) = J2 (cr ies ( r )hf-\r)+Cf Tes (r)hf )+ \r)) . 

i 

(4.26) 

Lastly we discuss how to construct the remaining fields 
using the gauge equations using the hierarchical scheme 
outlined in Table I. This is achieved by noting that the 
gauge equations (2.13)-(2.16) are for the retarded field. 
On the worldline we can write ftW = ZiW p . The 

remaining residual fields can be obtained by substituting 
this split into the gauge equations and rearranging for 
the h^ les . 


V. NUMERICAL IMPLEMENTATION AND 
RESULTS 

The self-force experienced by a particle moving along 
a fixed geodesic of the background Schwarzschild space- 
time was first calculated in the Lorenz gauge by Barack 
and Sago [30]. In calculating the retarded field they used 
a time-domain implementation for the l > 2 modes and 
used a frequency-domain method to calculate the mono¬ 
pole (Z = 0) and dipole (Z = 1) modes [32, 37]. They con¬ 
structed the self-force by projecting the tensor-harmonic 


modes of the retarded field onto a basis of scalar har¬ 
monics and regularizing using the standard mode-sum 
scheme. Lorenz gauge calculations were later extended to 
generic bound orbits in Schwarzschild spacetime [31, 33- 
35, 48, 49], 

In this section we detail how to compute, in the fre¬ 
quency domain, the Lorenz-gauge self-force along a cir¬ 
cular geodesic using the effective-source method we have 
developed above. Before giving the algorithm for the 
computation we briefly discuss how we construct numer¬ 
ical boundary conditions in order to solve for the retarded 
homogeneous metric perturbation. 


A. Numerical boundary conditions 


For the radiative modes (u> ^ 0) the asymptotic bound¬ 
ary conditions for the retarded field solutions are given by 
Eq. (2.17). In practice we cannot place the boundaries of 
our numerical domain at r* = ±oo. Instead we construct 
boundary conditions at a finite radius by expanding the 
asymptotic boundary conditions in an appropriate series. 
For the radiative modes we use the expansions: 


h {i) -(r in ) = e~ iUmr * 1 b k( r in - 2 M) k 
k =0 

fc+ax l 

fc (i)+ (r 0ut ) = J] 


k=0 


(5.1) 

(5.2) 


where r } n / out = r*(r in / out ). How the boundary locations, 
rin/out; and the truncation values, Zq[j ax , are selected in 
practice will be discussed in the algorithm section below. 
The series coefficients a k , b\. are found by substituting the 
above expansions into the field equations (2.10) and solv¬ 
ing for the resulting recursions relations. For brevity we 
do not repeat these relations here; they can be found in 
Appendix A of Ref. [33] . The recursion relations determ¬ 
ine the a l k>0 , b l k>0 in terms of the first coefficients a l 0 ,b‘ 0: 
respectively. By selecting appropriate linearly independ¬ 
ent vectors of these leading coefficients we construct a 
basis of linearly independent solutions that span the solu¬ 
tion space for the field equations. For example, for the 
odd radiative modes we have, once the gauge equations 
are employed, a solution space with 2 degrees of free¬ 
dom, i.e., we must solve for the h and h,( 10 \ For the 
outer homogeneous solutions the two basis are formed by 
setting {ao,aJ 0 } = {1,0} and {a^aj 0 } = {0,1}. Simil¬ 
arly we can repeat this with the {Z>o,Z>() 0 } f° r the inner 
solutions. 

In this work, although analytic solutions are now 
known [35], we opt to solve for the even static modes 
numerically as we already have code to do so. For these 
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modes the numerical boundary conditions take the form 
fe” ax 

h (i) -(n „)= bl(r in - 2M) k , (5.3) 

fc = fc~ in 
£.+ 

max 7 . —71 

S“> + (rou.)= E t gr °" (5.4) 

*=*+„ r ° ut 

How the truncation values fc± in are selected and the form 
of the recursion relations for the a l k , a k , b\ is again given 
in Ref. [33]. The log term in Eq. (5.4) is added to ensure 
the recursion relations have sufficient degrees of freedom 
to span the space of solutions to the field equations. 


B. Numerical algorithm 

The following steps describe how we calculate self-force 
in practice via our frequency-domain effective-source ap¬ 
proach. 

1. Choose a radial grid to store the values of vari¬ 
ous fields on. In general we require high resolution 
near the particle, and lower resolution far from the 
particle. Though our chosen window function is not 
formally compact, within our numerical procedure 
it is effectively compact. It is inside this effectively 
compact region that we need high resolution. In 
general we choose our window function to be effect¬ 
ively zero outside the region (ro — 2M, ro + 2M). In¬ 
side this region we find a grid spacing of M/ 10 suf¬ 
ficient (we pick the grid so that it includes r = ro). 
Outside the (effective) support of the window func¬ 
tion we use a grid spacing of 2 M. 

2. For radiative modes (m ^ 0) and even static modes 
(6 = even > 2,m = 0) construct numerical bound¬ 
ary conditions at r = 7' ou t and r = ri n using 
the recursion relations in Appendix A of Ref. [33] . 
For each £m- mode there will ben/ inhomogeneous 
fields to solve for - see Table I - and correspond¬ 
ingly nf sets of boundary conditions for the homo¬ 
geneous fields will be constructed as described in 
Sec. V A. 

3. For a given £m- mode solve for each basis of homo¬ 
geneous solutions and store the values of the fields 
and their derivatives on the preselected radial grid 
points. 

4. For the odd static (£ = odd, to = 0) modes and 
the monopole {l = m = 0) the values of the 
(in)homogeneous fields and their derivatives can be 
computed analytically. See Appendix D for an ex¬ 
plicit overview of the calculation for the monopole 
mode. 


5. For the given £m-mode compute the effective- 
source vector and store the results on the radial 
grid. 

6. At each grid point invert the matrix of homogen¬ 
eous solutions defined in Eq. (2.19) (formed from 
the previous stored results) and multiply it by the 
source vector to form the integrand of Eq. (4.25). 
Store the resulting values of the weighting coeffi¬ 
cient integrands at each point on the radial grid. 

7. Interpolate the weighting coefficient integrands us¬ 
ing standard cubic spline techniques. Numerically 
integrate the integrand as described by Eq. (4.25). 
The regular radial metric perturbation fields and 
their radial derivative are then computed via 
Eqs. (4.24) and (4.26), respectively. 

8. The gauge fields are then constructed as discussed 
at the end of Sec. IV B following the hierarchical 
structure given in Table I. 

9. The metric and its derivatives are constructed using 
the formulae given in Appendix A 6. The self-force 
is then constructed via Eq. (3.1). 

For comparison we also compute the retarded field with 
the method described in Sec. IIC. The first four steps of 
the algorithm in this case are the same as above. Then for 
the fifth step we use Eq. (2.22) to construct the retarded 
field weighting coefficients. This step only requires know¬ 
ledge of the homogeneous fields and the sources given in 
Appendix B. The retarded solutions are then constructed 
using Eq. (2.21). Finally we compute the self-force with 
both the standard mode-sum prescription described in 
Sec. Ill A, which relies on projecting the retarded tensor- 
modes onto a basis of spherical-harmonics before reg¬ 
ularization, and the tensor mode-sum prescription we 
present in Sec. VI. Note that only the radial compon¬ 
ent of the self-force requires regularization as it is the 
only component with non-zero regularization paramet¬ 
ers (see Eq. 3.12). Correspondingly, the contributions to 
both the t- and (^-components of the self-force converge 
exponentially in both i and l (the 0-conrponent is zero 
by symmetry). 


C. Results 

Using the above algorithm we can compute the resid¬ 
ual metric perturbation at r = ro- Using the residual 
field at the particle we can compute h^} u and we find our 
results agree with previously published results to relative 
accuracy of 10”'. Taking a radial derivative of the resid¬ 
ual metric perturbation at the particle we can compute 
the radial self-force without any further regularization re¬ 
quired. For the radial self-force we find agreement with 
the previous published results to a relative accuracy of 
10~ 6 - see Table III. In Fig. 2 we plot the residual field 
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for the = (2,2,1) field for a particle orbiting at 

ro = 6 M. 

A key feature of our procedure is that we only ever 
work with tensor-harmonic modes in constructing the 
self-force. This is in contrast to the standard mode-sum 
scheme whereby the tensor-harmonics of the retarded 
field are projected onto a basis of scalar harmonics before 
regularization. This projection, though straightforward, 
is cumbersome to implement (see, e.g., Appendix C in 
Ref. [30]). Furthermore the coupling between the tensor 
and scalar modes means that in practice to calculate £ max 
scalar modes one needs to calculate f? max = W + 3 
tensor modes. With our prescription this is not neces¬ 
sary. 

In Fig. 3 we show the convergence of the tensor mode- 
sum for the regular contributions to h^ u and F r . The 
punctures we use in this work are sufficiently regular that, 
for high-£ the contributions to h^ u and F r drop off as £~ 4 
and £~ 2 , respectively. 

In the next section we show how, by taking the limit 
to the worldline in our effective-source procedure, we can 
formulate a tensor-mode mode-sum scheme. 



ro /M 

this work 

Akcay et al. [18, 34] 

rel. diff. 

h K 

,L UU 

6 

-1.0471852(4) 

-1.0471854796(1) 

2 x 10” 7 

F r 

6 

2.4466487(8) x 10“ 2 

2.4466495(4) x 10“ 2 

3 x 10“ 7 

/)« 

,L UU 

10 

-0.48925802(2) 

-0.48925800172(4) 

4 x 10 f< 

pr 

10 

1.3389466(3) x lO" 2 

1.3389465(7) x 10“ 2 

3 x 10~ 8 


Table III. Sample results for orbits with r o = 6 M and 
ro = 10 M computed with £ max = 40. The relative differ¬ 
ence between our results and previously published data is 
small, being always less than 3 x 10 -7 . These results were 
computed using the higher-order punctures we provide online 
and by numerically fitting for the higher-order regularization 
parameters in order to speed up the convergence of the sum 
over tensor Amodes. Numbers in brackets denote the estim¬ 
ated error in the final digit of the corresponding result. Note 
that the results in this table have been adimensionalized, i.e., 
h^ u here = (M/^i)h^ u and F r here = (M/fi) 2 F r 


VI. MODE-SUM REGULARIZATION WITH 
TENSOR-HARMONIC MODES 



Figure 2. Sample results for the £ = 2, m = 2 mode for a 
particle orbiting at ro = 6 M. Shown are the h ^ and (scaled) 
h 1 ' 1 ' 1 metric perturbations for both the residual and retarded 
fields. At r = 6 M, the upper solid (red) curve shows the resid¬ 
ual field h^ Tes . The dashed (purple) curve shows the retarded 
field /i < ' 1)ret . Far from the particle the two coincide. Similarly, 
the lower solid (blue) curve shows /i ( ') res an d the dot-dashed 
(orange) curve shows h l - 7 ' >ret . The inset shows h^rl lea near the 
particle. With the punctures we present in the main text the 
residual fields are C 1 at the particle and correspondingly, as 
the inset shows, the second radial derivatives of the residual 
field are discontinuous. As the residual fields are C 1 at the 
particle the self-force can be directly computed from their de¬ 
rivatives. The punctures we provide online give C 2 residual 
fields at the particle which acts to improve the convergence 
rate of the mode-sum, as shown in Fig. 3. 
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In addition to their use in the effective source approach 
described here, the puncture fields may also be used 
to improve the efficiency of the traditional mode-sum 
scheme. In the standard mode-sum prescription, regu¬ 
larization is achieved through mode-sum formulae such 
as Eqs. (3.10) and (3.11), which take the form 


Figure 3. Convergence of the tensor Amode contributions 
to (the adimensionalized) h uu and F r for a particle orbiting 
at ro = 6 M. The punctures we present in this article in 
Appendix C are sufficiently regular that the contributions to 
huu and F r drop off as £~ . Online we given higher-order 
punctures that improve the rate of convergence for h U u to 
i~ 4 . 


TO _ 

b P ~ 


oo 

X[‘ 


e =o 


+ Dp, 

(6.1) 

recall, 

an £ 


sub/superscript denotes the scalar-harmonic multipole 


contribution (summed over m). The regularization para¬ 
meters Fl ^ and F ^ are analytically-derived functions 
of the instantaneous worldline. Provided the Detweiler- 
Whiting singular field is used in their derivation (and 
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the retarded field inodes F^ let are in Lorenz gauge), the 
D M term vanishes. Each term in the sum goes like £~ 2 
and so the partial sums converge as £^ ax . One can de¬ 
rive additional higher-order regularization parameters to 
accelerate this rate of convergence. For example, by sub- 

F [ 2 ] 

tracting an appropriate term of the form - - one 

finds that the terms in the sum now fall off as £~ 4 . One 
can continue in this way to higher orders, where the order 
£~ n term has the form 

OO 

_ 2£ ~+i _ = n ('6 2') 

(2£-n+l)(2£-n+3)...(2£+n-l)(2£+n+l) ’ v ' ' 

1=0 

for n even. The fact that the infinite sum over £ van¬ 
ishes is important as it guarantees that the subtraction 
of higher order regularization parameters does not affect 
the numerical result (other than accelerating the rate of 
convergence); equivalently, these higher order terms can 
be seen to come from pieces of the Detweiler-Whiting sin¬ 
gular field which vanish when evaluated on the worldline. 

In the scalar field case, this mode-sum formula is a 
natural choice as one can choose to work with scalar 
spherical harmonics labelled by £ when solving the field 
equations. In the gravitational case, the tensor spherical 
harmonics are a more natural choice and a numerical cal¬ 
culation typically produces tensor harmonic modes (la¬ 
belled by £) for the retarded field. Despite this fact, exist¬ 
ing calculations have relied on a scalar-harmonic mode- 
sum formula of the form given in Eq. (6.1). As a res¬ 
ult, a necessary step in the regularization procedure is 
the projection of the tensor-harmonic modes, E^ ret onto 

scalar harmonic modes F,( ret . This is undesirable for at 
least two reasons: (i) the projection involves cumbersome 
mode coupling formulas which have to be derived on a 
case-by-case basis, (ii) a given scalar-harmonic mode £ 
couples to several several tensor-harmonic modes (up to 
£ ± 2 for the metric and higher for some of its derivat¬ 
ives). This second point means that in order to obtain 
a given number of scalar £ modes, one actually has to 
compute several higher tensor-harmonic £ modes of the 
retarded field, and these are then lost during the projec¬ 
tion. Given that the cost of computing a given retarded 
field mode grows quadratically with £, this turns out to 
be quite a significant increase in computational cost. 

Fortunately, it turns out that the projection onto scalar 
harmonics is unnecessary; in this section we will derive 
tensor-harmonic regularization parameters which com¬ 
pletely eliminate the need for scalar harmonics. This 
addresses both issues mentioned above and produces a 
much simpler, more accurate and computationally effi¬ 
cient result. 

First, we consider what form a tensor-harmonic mode- 
sum scheme should take. The ^-dependence of the term 
of order £~ n will be given by 


for to > 0 and n integers. When m = 0 we can see 
that this reduces to the scalar-harmonic case, Eq. (6.2), 
as expected. For m > 0 and n > 2 the infinite sum 
over £ of any of these terms is zero, meaning we are free 
to add them without modifying the final result (other 
than accelerating convergence). In Eq. (6.2) this was 
only true when the sum starts at £ = 0, which would not 
be appropriate for tensor harmonics. For our generalised 
expression, Eq. (6.3), this holds for the sum starting at 
any value in the range 0 < £ < to. In practice we will 
have to equal to the value of m! used in the punctures 
and n will be determined by the power of p appearing in 
the singular field. 

Examining the puncture fields in Appendix C we can 
see they are already written in a form where this £- 
dependence is manifestly apparent. The task of produ¬ 
cing tensor-harmonic regularization parameters is there¬ 
fore merely a matter of reconstructing the £ modes of the 
singular metric perturbation using the expressions given 
in Appendix A 6, summing over to and then evaluating on 
the worldline. It is most convenient to do so in the (a, /3) 
coordinate system, i.e., using the punctures without the 
Wigner-D rotation matrices and evaluating at a = 0, as 
then only a small number of m' modes must be summed 
over. The only caveat is that this yields the components 
of the metric in the (a, ft) coordinates. We must there¬ 
fore also include a factor of the Jacobian from (a,/3) to 
(0, p) coordinates. This Jacobian is given by 


da 

1)0 

da 

dp 

dp 

dO 

dp 

dcj> 


— cos a sin p 
y/l — sin 2 a sin 2 fj 


0 , 


(6.4a) 


COS P ~ 1, 


(6.4b) 


— cos/3 

sin a\J 1 — sin 2 a sin 2 p 
— cos a sin /3 ^ 
sin a 


1 

sin a’ 


(6.4c) 


(6.4d) 


Note that because of the factor of =A_. appearing here, 
it is important to multiply by the Jacobian before taking 
the limit a — > 0. 

Since we are interested in computing the metric per¬ 
turbation and its derivative (for the self-force), we require 
mode-sum formulae for all components of the metric per¬ 
turbation and its derivative, i.e., 


h 


R 

fllS 




(6.5) 


and 


h 


R 


e te 


£=0 


(2£ + l 


( 6 . 6 ) 


A 


m,n 


2 n - 2m (2£+l)(£-m+l) 2m 
(2£ - 2m + n + 1) (£ - to + § + f ) 2m _„ 


(6.3) 


Here, the retarded-field £ modes are computed in the 
usual way from numerical data in the (0, p) coordinates, 













t 

KT= E ^( r o.f.°). ( 6J ) 

m——£ 

where the h?™ are constructed by combining the h ^ with 
the spherical harmonics; explicit expressions are given 
in Appendix A 6. The regularization parameters for the 
metric perturbation are computed using 


fc[°] = 


E 


P £m! 


dx M dx~ l 
dx^ dx u 


X—Xq 


( 6 . 8 ) 


where Xq denotes the point on the worldline, i.e. r = ro 
and a = 0 = /3. The only non-zero contributions in our 
case come from ml = 0 in the scalar sector [i = 1,3,6), 
m! = ±1 in the vector sector (i = 4,8), and to' = ±2 
in the tensor sector (i = 7,10). The regularization para¬ 
meters for the radial derivative of the metric perturbation 
are computed using 


hl $,T= E | 

m'=—£ 


dx M dx v ci z, 
dx* dx” u r ,l n'v 


Pirn' 


(6.9) 


where, again, the only non-zero contributions in this case 
come from m! = 0 in the scalar sector (i = 1,3,6), 
m! = ±1 in the vector sector (i = 4,8), and m! = ±2 in 
the tensor sector (i = 7,10). Finally, the regularization 
parameters for the ip derivative of the metric perturba¬ 
tion are computed using 


ft[°l 



dcx. ff ( dx^ dx u lP £m'\ 
dtp a y dx^ dx v 'V /y/ J 


X—Xq 


( 6 . 10 ) 


where the only non-zero contributions in this case come 
from m! = ±1 for hf r tm ' (i.e., i = 2) and m! = (0,±2) 
for hH m ' (i.e., i = 5,9). Evaluating these with the punc¬ 
tures given in Appendix C yields the following tensor- 
harmonic regularization parameters: 
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(6.11b) 

(6.11c) 

(6.lid) 

(6.lie) 
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h [ ~ 1] - T 
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(r 0 - M) 
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(6.Ilf) 
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32M(£ + 2/C) 


7r(r 0 -3M) 1 /2( ro _2M) 1 /2 
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/i 101 =- 

l0 tr,<p 
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; [ 0 ] _ 16[(r 0 - 2M)g - (r 0 - 3M)/C] 

^ 7r(r 0 — 2M) 1 / 2 (r 0 — 3M) 1 / 2 ( 1+ 

(6-llq) 

where, recall, Ai and A 2 are given by Eqs. (4.9) and 
(4.10), respectively, and where we have indicated with a 
subscript the cases (/i/^ 1 /, /4e r and h\p^} r ) where a terms 
is only non-zero above some minimum value of l. In all 
of the above equations, to simplify the presentation we 
have omitted an overall factor of the small mass /i. 

Finally, we note that these expressions can be com¬ 
bined to produce tensor-harmonic regularization para¬ 
meters for the redshift invariant h lll/ u ,1 u 1 ' and the radial 
component of the self-force. Doing so, we find 

H lo] = 


7 W r o+ L o 
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/C 


8/xM(6r 0 - 17M)/C 
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420p,M 2 /C 

nr 0 (r 0 — 3M) 3 / 2 (ro — 2M) 1 / 2 ’ 
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(6.14) 


Note that in giving these parameters we have rewritten 
Ai and A 2 in a form which highlights the fact that Jf[°l 
and .Fjgj both match their scalar-harmonic counterparts, 
Eq. (3.12), with the exception of higher order terms in 
1 /£. Since these terms vanish when summed from 1 = 0 
to infinity they have no impact on the final result and can 
be ignored in practice. Importantly, this is not the case 
for F[_ 1 j which differs from its scalar-harmonic version. 
The difference is in the presence of the second and third 
terms, and arises from the fact our mode sum expression, 
Eq. (6.9), starts at £ = 0 while it should start at i = 1 
for the vector sector and at l = 2 for the tensor sector. 
However, since this term has different limits on either side 
of the worldline, it vanishes upon averaging the left and 
right radial limits. As such, we see that in this case regu¬ 
larization can be achieved without projection onto scalar 
harmonics, by using the scalar-harmonic regularization 
parameters combined with an averaging procedure. 


VII. CONCLUDING REMARKS 

In this work we have developed a frequency-domain 
application of the effective source approach to computing 
the self-force on a point mass in a curved background 
spacetime. This new method builds on previous work 
which studied the case of a scalar-field toy model [28], 
extending it to the more physically relevant gravitational 
case. 

With a numerical implementation for the case of a 
circular orbit in Schwarzschild spacetime, our results 
demonstrate that the method can reliably produce accur¬ 
ate numerical results for the regularized metric perturb¬ 
ation with modest effort and computational cost. While 
this is not particularly important in a first order calcula¬ 
tion — the traditional mode-sum method, for example, 


can already produce comparable results with similar, or 
better, computational efficiency — the primary goal of 
our approach is to develop a set of methods which will be 
applicable to a second-order self-force calculation. Our 
results provide two key components in that regard: 

1. Our numerical scheme for solving the sourced field 
equations in the frequency domain carries over im¬ 
mediately to second order. The only change will be 
that the source will be a more complicated function 
involving the first-order metric perturbation. 

2. The source for the second order field equations is 
most efficiently written in terms of the first order 
Detweiler-Whiting regular field in an extended re¬ 
gion near the worldline. Such an approximation is 
exactly the output from our first-order calculation. 

In addition to addressing several important aspects of 
a second-order self-force calculation, the tensor-harmonic 
regularization parameters we derive in Sec. VI can also 
be used to improve the computational efficiency of a 
first-order mode-sum self-force calculation by avoiding 
the need for a cumbersome and wasteful projection 
onto scalar harmonics. It is interesting to note the 
close relation between the tensor harmonic regularization 
parameters and those one would obtain using a scalar- 
harmonic decomposition. In particular, provided one 
computes an average of either side of the (radial) limit 
to the worldline, we have found that scalar-harmonic 
regularization parameters may be used in place of their 
tensor-harmonic counterparts. We anticipate that this is 
more than merely a coincidence - in a future work we will 
investigate whether a similar result holds in more general 
cases. 

There are several future directions in which our res¬ 
ults may be extended. Most important is the application 
of our approach to the calculation of conservative effects 
from the second order gravitational self-force [9, 24, 25]. 
In addition to this, it may be interesting to study exten¬ 
sions of the approach beyond circular orbits, to the Kerr 
spacetime and to radiation and Regge-Wheeler gauges. 
With a view to identifying other important second-order 
effects, it may also be interesting to incorporate our 
method into an orbital evolution scheme which makes use 
of a two-timescale expansion of the equations of motion 
[50]. Such a scheme would likely provide a compelling 
balance of computational efficiency and faithfulness to 
the underlying physics of the EMRI problem. 
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Appendix A: Functions on the two-sphere 

When dealing with functions on the two-sphere, there 
are a wide number of possible conventions. Our conven¬ 
tions, which are consistent with those of Mathematica 
[51] are summarized in this Appendix. 


1. Scalar spherical harmonics 


The associated Legendre polynomials may be defined 
in terms of derivatives of the standard Legendre polyno¬ 
mials, 


pr(x) = (-ir a-* 2 ) 


_ ™2\m/2 


dx ? 


-Pt{x) [m > 0], 




(Ala) 

(Alb) 


where we have included the Condon-Shortley phase 
factor (—l) m and where the Legendre polynomials satisfy 
the Legendre equation 


(I-* 2 ) 


2 d 2 P t {x) dP t (x) 


— 2x- 


+ £{£ + l)P e (x)=0. (A2) 


dx 2 dx 

We now define the scalar spherical harmonics as 


Y (rn (e,ip) = 


{21 + 1) {£ — m)\ 
47t {£-\-m)\ 


P™(cos0) e im,p , (A3) 


where —£< m < £. Note that £ and m are merely 
labels and we will raise and lower their position freely 
to wherever they get in the way the least. The scalar 
spherical harmonics are orthonormal, 

/> 27T /* 7T 

/ / Y lm {9,ip)Yt, ml (e,<p)d£l = 6 U '& m m', (A4) 

Jo Jo 

where dfl = sin ddddip. They also satisfy 

y/ m (0, <p) = (-1 rw,_ m (0, <P), (A5) 

and the completeness relation 

OO £ 

Y Y Y tm{0,v)Y!J m {e',y') = 8{cos8-cosd')5{p-ip'). 

£—0 m——£ 

(A6) 


Since the scalar harmonics form an orthonormal basis, 
an arbitrary scalar function can be expanded in spherical 
harmonics, 


oo £ 

f(8><P)=Y E f e m Y tm{0, <fi), (A7) 

£—0 m——£ 

where the coefficients are given by 

/*27T /»7T 

fim= / f{6,tp)Yl m (6,tp)dn. (A8) 

Jo Jo 

At the pole, 0 = 0, only the m = 0 spherical harmonics 
are non-zero and (A7) becomes 

f{®,v) = Y\j 2 ^tr fm ' (A9) 

<=o 


where 

lop _j_ i r 27r r n 

feo = y 47r J J f{a,/3) P e {cosa)dn. (A10) 
Similarly, the 8 derivative only requires modes m — —1,1: 


{def){ 0,<p) 


E 

t—0 


£(£+l)(2£+l)^_ lv 


167r 


{e~ lv fc,-i 


the second 6 derivative requires m = 0, ±2: 


e iv fn), 

(All) 


(d sg f){0,<p) 


OO 


E 


( 2 £ + 1 ) 

647T 


2£{£+l )fto 


+V(* - !)£{£ +1){£+ l){e~ 2i Pfe -2 + e 2 ^f e2 ) , 

(A12) 


and so on; for n derivatives with respect to 8 we need 
modes m = —n, —n + 2, • • • , n — 2, n. 


2. Vector spherical harmonics 

The vector spherical harmonics fall into two categories, 
those of even parity (£ + m even) and those of odd parity 
(£+rn odd). These categories reflect a difference in beha¬ 
viour under the parity operation (0, ip) -A (7r — 0, ip + 7r); 
the even parity harmonics are invariant under this trans¬ 
formation while the odd parity harmonics change sign. 

The even-parity vector harmonics are defined by 

= [£{£ + 1 ))-^ 2 D A Y lm , (A13) 

and the odd-parity harmonics are defined by 

X^ 1 = -[£{£ + 1 )]-^ 2 e A B D B Y em , (A14) 

where Da is the covariant derivative and cab is the 
Levi-Civita tensor associated with the metric Dab = 
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diag(l, sin 2 0 ) on the two-sphere (i.e. eg v = sin 0 , e v g = 
— sin0, egg = 0 = e vv ). Explicitly, the components of 
the vector harmonics are 


Z e e m = [£(£+1 )\~ l / 2 dgY tm , 

X% m = -[1(1+ \)]-^-^d v Y tm , 

Z t ™ = [l{i + l)]- 1 l 2 d v Y tm , 

X*™ = [£(£ + l )]" 1 / 2 sin 0 5 gY lm . (A15) 

The vector harmonics satisfy the orthonormality rela¬ 
tions 



f A^ m (0, ip) Xp m , (0, if) dfi = 5w Sr, 
JO 

r ¥>) Ztm' ( 0 , <p) dn = 8u> S„ 

J 0 

r X l r{0,<p)Z$* m ,[p, ip) d£l = 0. 


Jo Jo 

They also satisfy 


(Al 6 a) 

(Al 6 b) 

(A16c) 


XTifiM = (-1 ) m X%- m (0,ip), (Al7a) 

ip) = (-1 ) m 4’- m (0, ip). (Al7b) 


These definitions are consistent with [52] and with [53] 
(apart from the inclusion of the prefactor [£{£ + 1)]~ 1//2 
which ensures orthonormality) and relate to those of [54] 
through the conversion Z 1 / 1 -» [£(£ -fl)] 1 / 2 y| m , X^ -4 
[£{£+! )] l /*X l j?. 


3. Tensor spherical harmonics 

The tensor spherical harmonics again fall into two cat¬ 
egories, those of even parity (£+m even) and those of odd 
parity (£+m odd). The even-parity tensor harmonics are 
defined by 


+ \£{£ + 1) sin 2 0 Y em , (A20c) 


X tm 

06 ~ ~ 


X £m _ 

6<p — ~ 


\r£m 



. n cot 0 du> 

sin 0 L 


Y Em , (A20d) 


1 


SJ<pif j sin 


9 (*- 2 )i' 

2 

sin 0 

(£ + 2 )! 



2 sin 0 . 

+ sin 0 cos 0 dg 

dg v - cot 0 d v 


Y 


£m 


(A20e) 


d e Y tm . (A20f) 


The tensor harmonics satisfy the orthonormality rela¬ 
tions 


/*27T /* 7T 

/ X%{6,ip) X AB t (0, ^) cM = 8w 

Jo Jo 

/*27T /»7T 

/ / z%(e,<p)z^t(d 1 < p )dn = 6u' 

Jo Jo 

/* 27T /»7T 

/ / x5s(0, VJ )z^f(0 )V )dn = o, 

JO JO 

and the identity 

n AB z% = o = n AB x%. 


Smm' * (A 21 a) 
<W,(A21b) 
(A21c) 


(A22) 


They also satisfy 

*aj? (*,¥>) = (-l) m ^- m (0,^), (A23a) 

= (A23b) 


These definitions are consistent with Thorne [52] and re¬ 
late to those of [53] through an orthonormality factor, 


r/Urn 
AAB 


o (1—2)! ^ xr£m ytm 
Z (t+2)! 1 ABi a AB 


-A 


9 (1-2)1 
Z (t+2)! 



X 


AB- 


4. Rotations 


Z tm _ 

AB — 


'o (l- 2 )! ~ 

. (€ + 2 )!_ 


[D a D b + *£(£+1 )D ab ]F^, 

(A18) 


and the odd-parity harmonics are defined by 


X £m _ 

AB ~ ~ 


o (l— 2 )! ~ 

. (€ + 2 )! 


e {A c D B) D c Y tm . 


(A19) 


Explicitly, the components of the tensor harmonics are 


1 


II 

9 (^ — 2 )! 
(£ + 2 )! 

2 

5 0e + ^(f + i)]y^ m 

'£m 

"0 — 2 )i" 

1 

2 

5^-cot 05^™, 

6(p 

“(* + 2 )! 


II 

£ $- 

" (€- 2 )!' 
(€ + 2 )! 

1 

2 

d vv + sin 0 cos 0 5# 


(A20a) 

(A20b) 


Under a rotation of the coordinate system which is 
represented by the Euler angles a, /3, 7 , the spherical har¬ 
monics components transform according to 

t 

fem{0,<p)= ^2 D iim'( a ^,l)fem'(0 , ,p'), (A24) 

m'=—£ 

where 3, 7 ) is the Wigner-D matrix [55]. Here, 

we use the convention that the Euler angles correspond 
to a z — y — z counterclockwise rotation and our con¬ 
vention 2 for D t rnm ,{aj 3 , 7 ) is consistent with Rose [56]. 
Using these conventions, the Wigner-D matrix satisfies 

D e mi m 2 (<Yf 3,7) - e-™i™^Di im 2 (0,/3,0). (A25) 


2 This convention is different from that of Mathematica [51] and 
Wigner [55]. Our , (ol,/3 , 7 ) is related to theirs by a change 
in the signs of m and m' [56]. 
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The vector and tensor harmonics also transform in a sim¬ 
ilar way [57], i.e., 

xT(o,v) = 

r)r A ' 1 

E 7 )^( 0 ',^), (A26a) 

m'——£ 

zT{o,v) = 

a t a ' t 

-g^A E (A26b) 

m'——i 

x%(fiM = 

f)r A ' f) t b ' JL , 

-Q^A^B E (A26c) 

m'=—£ 

z%(e,v) = 

f)r A ' f) t b ' , 

E nLm'^^z^ie',?'), (A26d) 

m'——£ 

which is equivalent to stating that the vector and tensor 
harmonic components of a tensor transform according to 
Eq. (A24). Finally, we note that the Wigner-D matrix 
relates to the spin-weighted spherical harmonics: 

DL(*,P, 7) = (A27) 

which for the spin-0 case gives a relation to the scalar 
harmonics, 

Dl l 0 (a,AO) = y^y; ra (Aa). (A28) 


5. Tensor harmonic basis in Schwarzschild 
spacetime 

Barack and Lousto [36] use the above bases of scalar, 
vector and tensor harmonics to construct a basis of har¬ 
monics for the components of a symmetric rank-2 tensor 

defined on a Schwarzschild background spacetime. 
This basis was later modified slightly by Barack and Sago 
[30] to improve the behaviour of some components near 
the horizon. In particular, they choose a basis of 10 fields 
in t — r space defined by 

p2~K p 7T 

42 = / r (*tt + f 2 trr)Y e * m dfi, (A29a) 

Jo Jo 

p2lV p 7T 

42 =/ / 2rft tr Y; m cm, (A29b) 

Jo Jo 

p2n p 7T 

42 =/ / rf (t tt - f 2 t rr )Y; m dn, (A29c) 

Jo Jo 

p2,7T p 7T 

42 = / / 2 [€(* + 1)] dfi, (A29d) 

Jo Jo 

p2n p 7r 

42 = / / 2 [*(* + 1)] 1/2 ft rA z£* n dn, (A29e) 

^0 ./0 


f ( 6 ) _ 


t (?) - 


t»27T /* 7T 


/o Jo r 

p2ir p 7T -i 


t AB n AB Y e * m dn, 


I o Jo 


, (l — 2 )! 

r n* + 2)!. 


1/2 


tAB ( 


ryAB* (~)AB(~) ryCD*\ 
^im ~ ManAt. 1 


CD^£ m j dn, 

(A29f) 

p'A'K pTT 

/ / 2[^+l)] 1/2 t 4A X/ ro *dO, (A29g) 

J 0 </0 

2 [f (£ + 1)] 1/2 t rA Xf^ dfl, (A29h) 


(8) 

(9) 

im 


tO JO 

p2n pTT 


to JO 


Aw) _ 

L im ~ 


r r - 


Jo Jo r 

[ (^ + 2)! J 


t AB (x AB * - Ui An n C DX^) dn, 

(A29i) 


where / = (1 — 2 M/r). The harmonics i = 1,... ,7 are 
of even parity, while the harmonics i = 8,9,10 are of odd 
parity. 

Barack and Sago represent this basis in terms of a set 
of 10 tensors defined by 


II 

^W+rXsi)y tm , 

(A30) 

y( 2 ) — 
^1/ 

^(s^ + s^Dy^, 

(A31) 

y(3) — 

jsW-rXWY*™, 

(A32) 

II 

+zfrsi), 

(A33) 

^ Ot 

II 

j^{5 r x m + zfrsi), 

(A34) 

y(6) — 

/J.U 

^n AB 5 A s B Y em , 

(A35) 

II 

r\z^ - \z A A n^ v ), 

(A36) 

ii 

00 s ^ 

-^x^+xfrsi), 

(A37) 

II 

oT ^ 

-771 (5;xfr+xfr5:), 

(A38) 

y(io) _ 

r 2 {X e ™ - i*W)- 

(A39) 


With the exception of * = 3, this basis is an orthonormal 
set in the sense that 


p 2 tt p 7T 

/ / V T 'V% ( i )/m Y^/ m ‘* dn = Sij 6w 5mm ', 

Jo Jo 

(A40) 

where rf K = diag(l, / 2 , r -2 , r -2 sin -2 0). For i = 3, the 
set is also orthogonal, but Y$ has a norm of / 2 . 

Finally, we note that Barack and Sago factor the coef¬ 
ficients 




y/2 


1, i = 1,2,3,6, 

(^Tl))- 1 / 2 , z = 4,5,8,9, 

((£-l)£(£+l)(£ + 2))- 1 / 2 , * = 7,10, 

(A41) 


out from the tensor harmonic fields, ftW, in order to 
make some of their expressions for, e.g., the field equa- 
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tions more compact. We likewise use these coefficients in 
Eqs. (2.7) and (4.5). 


6. Metric reconstruction 

Rebuilding the original metric perturbation, h , from 

the fields is straightforward. The necessary equation 
can be derived using Eq. (A29) above along with the fact 
that a trace reversal, h^ = h^ — ^g^vh, is equivalent to 
the interchange <t=> f 2 h^\ This gives 


OO l 

^ = ^EE ( A42 ) 

i m——i | 


where 

^ = + (A43) 

hir = m-^Y* 171 , (A44) 

<7 = f(r)~ 2 ( h - fhfi) Y lm , (A45) 

h tA = {hfrlzT - hf^xT) , (A46) 

= f{r) J f(e+1) 0^” - S '"W) , (A47) 

h% = r 2 n AB h^Y^ 

+r2 /M(^ (Z - " ^ C c^ab) 

-^L 0) (4b - ^cOu*)), (A48) 

and where the sum over £ begins at £ = 0 for the scalar 
sector (i.e., hj m , hj m , and /if m ),_at ^ = 1 for the 
vector sector (i.e., hf m , h\ m , /if m and /if m ), and at £ = 2 
for the tensor sector (i.e., /ij m and hj^). 


Appendix B: Field equations and retarded field sources 

The coupling terms in the frequency-domain field equation (2.10) are given by 

M^ U) h^ = ± (l - (ftW - h (5) - A (3) ) - £ (l - 


Mw u) h& = \ffhw + \r 

_ff fz 

2r 

= - A 


»-5< 2 >) + y r 2 .>-s<« 

/j(5) _ //j(3) _ 2/h*- 6 ^ , 

/j(l) _ /j(5) _ 4 _ ^(3) + /j(6)^ 


/E 

2r 2 


(4 2 ) - fcW) 


M^ U) h^ = \r 

4 r 


VjJ 


- h^A +hA ] - hf 1 - 4(£+ l) 4 ^ (2) 

/ ’ * ’ *J 2 r z 


(3/r (4) + 2 h {5) - h (7) + £{£ + l)ft, (6) ) , 


A* ( V“ = £ [(l - “) ^ (i (1 > - A<’>) + ) 




m ( 6) w P = - 4 


2r 2 


/z. (1) - /i (5) - ( 1 - 


4M 


(/i (3) + /i (6) ) 


MV [i) hU=-±(hV + M*>) 


MW u) hW = -/' 


+ JiW - /$>] - ( 3 ~ h(8) + 2 h {9) - ^ (10) ) 


^ ( 10 ) y ) ^ ) =-^(fc (10) + A?i (9) ). 


l)/i (6) - 


(Bl) 

(B2) 

(B3) 

(B4) 

(B5) 

(B6) 

(B7) 

(B8) 

(B9) 

(BIO) 


where A = (£ — 1)(1? + 2). 
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The sources to the held equation (2.10) take the form 


where 


Jir) - -™aE a <f> s(r _ rn) / yrM ^ *=!>■ • ■ - 7 

’ f 2 a 0[r r °M F^ m *(7r/2,n v t), * = 8,9,10 ’ 


« (1) = fo/ r o, « (2) = 0 , a (3) = f 0 /r 0 , a (4) = 2 if 0 mQ. v , a (5) = 0 , 

a (6) = r 0 fl 2 , a (7) = r 0 ff 2 [t{t+ 1 ) - 2 m 2 ] , a (8) = 2 f 0 tt v , a (9) = 0 , a (10) = 2 *mr 0 fl 9 . 


(BH) 


(B12) 


Appendix C: Puncture functions for circular orbits in Lorenz gauge 


In this section, we give our explicit expressions for the Lorenz-gauge puncture fields, , for the case of a circular 
geodesic orbit in Schwarzschild spacetime. These punctures contain all pieces of the Detweiler-Whiting singular held 
necessary to compute the regularized components of the metric and its hrst derivatives. Written as tensor-harmonic 
modes in the (9, ip) coordinate system, the punctures are given by 


fWP _ r _ 
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+Ar 
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(27 — 1)(27 + 3) 


3tt M(r 0 - 2M) 1 / 2 (r 0 - 3M) 1 / 2 


+(2£ + l)|Ar 
+Ar ( 


Mr. 


1/2 


(r 0 - 2M)(r 0 - 3 M) 1 / 2 
/ 256[2(ro - 4M)(r 0 - 2 M)£ - (r 0 - 3M)(2r 0 - 7M)K] 


A 9 


37 rAf(r 0 - 2M) 1 / 2 ( ro _ 3M) 1 / 2 
1 80[(16r 2 - 96Mr 0 + 131M 2 )£ - 2(8- 52 Mr 0 + 81 M 2 )K] \ 


(2£ — 1)(27 + 3) 


3tt M(r 0 - 2M) 1 / 2 (r 0 - 3M) 1 / 2 


(CIO) 


where -D^ m / = U^ m /(7r, 7r/2,7r/2) is the Wigner-D matrix corresponding to the rotation from (a, (3) coordinates to 
(9,<p) coordinates and, recall, Ai and A 2 are defined in Eqs. (4.9) and (4.10), respectively. 
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Appendix D: Monopole contribution to the 
Lorenz-gauge metric perturbation 


In this section we calculate the monopole {t = 0) con¬ 
tribution to the retarded and residual metric perturba¬ 
tion at the particle. There is some subtlety to Lorenz- 
gauge monopole perturbations which we will highlight, 
but we refer the reader to the given references for more 
detailed information. We begin by providing the field 
equations and basis of homogeneous solutions for the 
monopole perturbation. 


1. Field equations and basis of homogeneous 
solutions 


The generic form of the held equations is given by 
Eq. (2.10). For the monopole perturbation, which has 
l = 0,m = 0 and ui m = 0, only the (i) = 1,3,6 modes 
are excited. The held equations can be further simpli¬ 
fied using the gauge equation (2.14) to decouple the h^ 
held — see discussion around Table I. The remaining held 
equations for are given by 


h (1) 

,rr 


(r - 4 M)/iW _ h w - f(rhW - /i< 3 >) J , 

(Dl) 



rhf - h (3) + ((4 M - r)/#> 



In order to display the basis of homogeneous solutions, 
let us define 


H = {M/n) {h tt ,h rr ,r 2 hgg = (rsinfl) 2 h vv ] 

(D3) 

where the second line derives from the metric recon¬ 
struction formulae (A43), (A45) and (A48), noting that 
loo = • The inverse relations are 


h (1 ) = 2 y/nn 1 r{h tt + f 2 h rr ), 

(D4) 

h (3) = Aypn : ir 1 r~ l h gg , 

(D5) 

h (6) = 2v/7r/r -1 y(/j tt - f 2 h rr ). 

(D6) 


A complete basis of homogeneous solutions to the two 
coupled monopole held equations (Dl), (D2) is given by 


[34] 

H A ={-f,f-\ 1} 


(D7) 

(D8) 

(D9) 


_ f M 4 M 3 f~ 2 (3M — 2r) M 3 1 
° l r 4 ’ r 4 r 3 J ’ 

Hd = {^ i w ( r ) + rP ( r )/ ln /- 8M " ln m) > 

f — [K(r)-rQ(r)f\nf-8M 3 (2r-3M)\n^] , 
^ [3r 3 - W(r) - rP(r)f ln / + 8 M 3 ln £] J 


(DIO) 


where 


P(r) = r 2 + 2rM+ 4M 2 , (Dll) 

Q(r) = r 3 - r 2 M - 2rM 2 + 12M 3 , (D12) 

W(r) = 3 r 3 - r 2 M - ArM 2 - 28M 3 /3, (D13) 

K(r ) = r 3 M - 5r 2 M 2 - 20rM 3 /3 + 28M 4 . (D14) 

By substitution, it is straightforward to verify that the 
set {Ha, H b , Hq, H D } are solutions to the homogeneous 
held equations (Dl) and (D2). 

When constructing the inhomogeneous monopole solu¬ 
tion it is important to ensure that it represents a particle 
with the correct mass-energy. That is, by Birkhoff’s the¬ 
orem, for the spherically symmetric monopole perturb¬ 
ation the solution must have the same geometry as the 
Schwarzschild solution with mass M for r < tq. For 
r > r o the solution must again be that of Schwarzschild 
geometry but with mass M + /iEq. As we now briefly 
discuss, perhaps the most natural method for construct¬ 
ing the inhomogeneous monopole perturbation does not 
satisfy this condition. 

First we note that the solutions Ha and H B are regular 
at the event horizon but approach nonzero constants as 
r —> oo. Conversely, He and H B are regular at infinity, 
but are singular on the horizon. It is therefore tempt¬ 
ing to construct the inhomogeneous ‘internal’ solution for 
r < ?'o as a weighted sum of the {Ha,H b } basis func¬ 
tions and an ‘external’ solution for r > ro as a weighted 
sum of the {Hc,H B } basis functions. This turns out to 
not give a solution that has the correct mass-energy [26] . 
Instead a correction term, AIAh, must be added [34] so 
that for the retarded field we have 


m o = Ai?ih 


CaHa + C b H b , r < r 0 , 
CcHc+CdHd, r>r 0 , 


(D15) 


where Ca, C b ,Cc, C b are constant weighting coefficients 
and 


AH ih = -C A (H A -H B ). (D16) 

A curious, but well understood, feature of the Lorenz- 
gauge monopole perturbation is that the tt component 
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approaches a nonzero constant at infinity [13, 37]. When 
computing gauge invariant quantities, care must be taken 
to account for this minor peculiarity of the Lorenz gauge 
as is discussed in Refs. [26, 58, 59] (for circular orbits) 
and Refs. [31, 60, 61] (for eccentric orbits). 

Before we proceed it will useful to define a matrix of 
homogeneous solutions by 


$ = 


h {1) 

n A 

-h {1) 

n B 


h (1) 

n D 

h {3) 

n A 

-h {3) 

n B 

h? 

h {3) 

n D 

/i (1) 
il A r 

-h {1) 

n B,r 

n ‘C,r 

h (1) 

n D,r 

ll A,r 

-h i3) 

u B,r 

h (3) 
n ‘C,r 

h {3) 

ll D,r 


(D17) 


where the matrix elements are constructed by applying 
the relations (D4) and (D5) to the basis of homogeneous 
solutions (D7)-(D10). 


2. Retarded solution for the monopole mode 

The retarded field at the particle is constructed via 
Eq. (2.21) where the sources are given by 


C c = 


2 [8Mr 0 - 3 rl - 12 M 2 + 24M(3M - r 0 ) In a] 
9My/r 0 (r 0 - 3 M) 


(D22) 



(D23) 


The monopole contribution to the retarded metric per¬ 
turbation everywhere in the spacetime is then given by 
Eq. (D15). 


3. Residual solution for the monopole mode 

The residual metric perturbation is constructed using 
Eq. (4.24) and Eq. (4.25). The effective sources, 
that appear in the latter equation are constructed by 
making the replacements JY 1 / 3 ) —i Wh^ 1 ^ 3 ^ in Eqs. (D4) 
and (D5) for and S^ 3 ) eff , respectively. The two ne¬ 

cessary punctures are given by Eqs. (Cl) and (C3) and 
W is the window function. The resulting effective sources 
are rather cumbersome so we will not display them but 
they are easily constructed using computer algebra pack¬ 
ages. 

Following Eq. (4.24) the (r-dependent) weighting coef¬ 
ficients can be solved for via 


J{1) = _8 Jl3) = JW 
ro Jo 

Following Eq. (2.22) the (r-independent) coefficients are 
given by 

(C A C b C c C d ) t = $- 1 .(0 0 jK 3 4 )f . ( Dig ) 


rb 

{CT Cb s Q? s C% S ) T = / ^(OO S (1)eff S (3)eB ) T dr, 

J a 

(D24) 

where for C A and Cb the limits on the integral are a = 
r,b = oo and for Cc and Cd the limits are a = 2 M, b = r. 
The monopole contribution to the residual metric per¬ 
turbation everywhere in the spacetime is then given by 


Explicitly the weighting coefficients take the form 


c 

\/ro(ro-3M)’ 

(D20) 

c 8M + (6M-2r 0 )ln/ 

3^/r 0 (r 0 -3M) 

(D21) 


Co(r) =A H ih + CT(r)H A (r) + C B (rY es H B (r) 

+ C c (ry es H c (r) + C£ s (r)H D (r). (D25) 

Here, although we are calculating the residual field, the 
AHih is still that given by Eq. (D16) with constant C A 
and Cb calculated via Eq. (D19). 
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